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Starting  Date:  09/10/2012,  End  Date:  06/09/2013  (ARO,  STIR  program) 

Chihoon  Lee 
Department  of  Statistics 
Colorado  State  University 

1  Statement  of  the  Problem  and  Summary  of  the  Results 

In  the  first  part  of  this  project,  we  study  stationary  solutions  of  constrained  stochastic  dif¬ 
ferential  equations  (SDEs)  driven  by  fractional  Brownian  motions  (FBMs).  Motivations  for 
this  study  originate  from  the  fact  that  such  constrained  processes  serve  as  an  approximation 
for  a  large  class  of  high-speed  communications  stochastic  networks  in  heavy  traffic,  such  as 
local  area  networks  and  internet  systems.  To  study  ergodic  properties  of  a  dynamical  sys¬ 
tem  with  memory,  we  investigate  the  ideas  based  on  Markovian  random  dynamical  systems 
and  establish  a  key  tightness  (in  time)  result  of  the  constrained  stochastic  processes.  In  a 
framework  of  Stochastic  Dynamical  Systems  (i.e.  infinite  dimensional  state  space  setting 
that  pertains  to  noise  process  with  memory),  such  a  tightness  result  essentially  establishes 
the  existence  of  the  stationary  solutions. 

In  the  second  part  of  this  project,  we  address  several  parameter  estimation  problems  for  a 
class  of  one-dimensional  reflected  SDEs  driven  by  FBMs.  We  note  that  parameter  estimation 
problems  are  usually  quite  challenging  in  the  physical  network  models  (with  or  without  long 
memory),  whereas  the  limit  SDE  models  can  be  much  more  tractable  for  statistical  analysis. 
We  derive  some  explicit  form  of  maximum  likelihood  and  method  of  moment  estimators  for 
the  drift  parameters,  then  obtain  their  desirable  statistical  properties,  which  are  based  on 
continuous  time  and  discrete  time  observations,  respectively. 

In  the  third  part  of  this  project,  we  justify  the  use  of  constrained  SDEs  as  approximation 
models  for  a  very  general  class  of  stochastic  networks  in  heavy  traffic  with  non-Markovian 
state  dependence;  long  range  dependence  can  be  considered  as  a  special  case  of  such  broader 
assumptions.  We  establish  heavy  traffic  limit  theorems  for  queue-length  processes  in  critically 
loaded  single  class  queueing  networks  with  state  dependent  arrival  and  service  rates.  The 
limit  stochastic  process  is  a  continuous-path  reflected  process  on  the  nonnegative  orthant. 
The  constrained  stochastic  differential  equation  driven  by  FBM  is  a  special  example  of  such 
limit  stochastic  processes.  We  also  give  an  application  to  the  classical  generalized  Jackson 
networks  with  state-dependent  rates. 

Overall,  the  conducted  research  would  advance  the  existing  knowledge  in  the  steady-state 
behavior  and  control  issues  of  stochastic  networks  operating  under  heavy  traffic  regime.  The 
work  is  expected  to  have  impacts  beyond  its  theoretical  contributions  in  applied  probability, 
with  advances  likely  to  be  influential  in  operations  research  and  defense-related  strategic 

1  Research  supported  by  US  Army  Research  Office  W911NF-12-1-0494. 
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Figure  1:  An  illustration  of  a  typical  stochastic  network  (at  left)  and  its  approximating  limit 
model  (at  right).  The  coordinate  Wi ’s  represent  the  weak  limit  of  the  suitably  scaled  queue- 
length  (or  workload)  processes  for  service  stations  of  the  underlying  network  system.  The 
approximations  give  parsimonious  mathematical  models  that  are  much  more  tractable  for 
statistical  inference  as  well  as  long  time  asymptotic  analysis. 

engineering.  Design  of  efficient  communication  networks  increasingly  relies  on  theoretical 
approximations  and  statistical  parameter  estimations  such  as  those  to  be  developed  here. 
While  the  research  is  motivated  by  stochastic  network  data,  the  methodologies  could  be 
utilized  in  applications  well  beyond  the  network  setting  and  likely  to  be  extended  and  adapted 
by  researchers  from  many  other  areas  of  science  and  engineering  where  long  range  dependence 
and  self  similarity  are  key  features. 

2  Stationary  Solutions  to  Constrained  SDEs  Driven  by  FBMs 

The  SDE  model  under  consideration  is  given  by 

W(t)  =  W( 0)  +  [  b(W(s))ds  +  BH{t )  +  rj(t )  e  X,  t>  0, 

Jo 

where  is  a  standard  fractional  Brownian  motion  (with  appropriate  dimension)  and  X  is 
a  given  (polyhedral)  domain.  Roughly  speaking,  the  nonnegative  process  rj ( • )  describes  the 
minimal,  continuous  correctional  term  with  finite  variation,  and  increasing  only  when  (kF(f)} 
is  on  the  boundary  dX.  In  other  words,  the  processes  {W(t),rj(t)}  solve  the  Skorokhod 
problem 

wit)  =  r  (V(o)  +  jf  b(w(s))ds  +  it)  ex,  t  >  o, 

where  T(-)  is  the  Skorohod  (reflection)  map.  Based  upon  such  a  representation,  we  say  that 
{ W(t)}t>o  is  a  stationary  solution  to  the  above  equation  if 

(w{t)-W{  0)-  f*  b(W  (s))ds  —  r)(t)\  ={BH{t))t>  o, 

\  Jo  /  t> o 
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where  the  equality  is  the  equality  of  all  finite  dimensional  distributions,  and  their  joint 
probability  distribution  does  not  change  when  shifted  in  time.  With  this  setup,  our  main 
goal  is  to  establish  the  existence  and  uniqueness  of  the  stationary  solution  under  the  suitable 
stability  condition  on  the  drift  term  b(-).  The  main  idea  is  based  on  studying  sequences  of 
empirical  measure  of  Euler  schemes  described  below. 

Fix  a  constant  step  size  7  >  0,  then  denote  (1147)fc>o  the  Euler  scheme  defined  by 
W0  =  x  G  X,  W(k+ 1)7  =  f  (Wk'y,  lb{Wkl)  +  Afc+1)  , 

where 

Afc+i  =  BH((k  +  1)7)  -  BH(k 7) 

is  a  fractional  Gaussian  noise  and  r(w,  •)  is  the  Skorokhod  map  T  evaluated  at  time  1  with 
initial  condition  w  G  X.  Then,  a  sequence  of  random  probability  measures  da))n>  1 

is  defined  on  the  Skorokhod  space  D(]R+,X)  by 

1  n 

P^\u,da)  =  - 

L  k=  1 


where  5  denotes  the  Dirac  measure.  This  random  measure  constitutes  our  basic  sequence  of 
approximations  for  the  stationary  solution  to  the  constrained  SDE  models.  Under  suitable 
stability  conditions  on  the  drift  vector  held  b(-)  and  regularity  of  the  Skorokhod  map,  we 
are  to  establish  that  (u ,  da)  for  large  n  and  small  7  are  natural  approximations  of  the 
stationary  solutions. 


A  fundamental  step  is  the  tightness  of  da))n>i. 

establish  that 


1 

sup  — 

n>  1  ^ 


|fU(fc_i)7|2  <  +00 

k= 1 


a.s. 


For  the  tightness  result,  we 


for  sufficiently  small  step  size  7  G  (0,  00).  The  stationary  and  ergoclic  properties  of  the 
fractional  Gaussian  noise  sequence  (A*,)fc>i  along  with  suitable  stability  properties  of  the 
Skorokhod  map  play  an  important  role  in  the  proof  of  the  above  estimate. 

Next,  we  identify  that  every  weak  limit  of  da))n>i,  denoted  by  z/^00,7^)(u;,  da ),  is 

indeed  a.s.  a  stationary  Euler  scheme  with  step  7  of  SDE  model.  More  precisely,  we  denote 
by  {W^'y'>}t> 0  the  stepwise  constant  continuous-time  Euler  scheme  defined  by  U77'1  =  !147 
for  all  t  G  [k'y,  ( k  +  1)7).  Then,  we  show  that  a.s.  i/(°°,7)(u;,  da)  is  the  distribution  of  an 
RCLL  process  =  (wS°°’7))t>0,  which  is  a  discretely  stationary  solution. 

Then  the  last  important  step  is  to  prove  that,  a.s.,  (z/^°°’7^(a;,  da))7>o  is  tight  for  the 
weak  topology  induced  by  the  topology  of  the  uniform  convergence  on  D(M+,  X)  and  that 
its  weak  limits  when  7  — >  0  are  stationary  solutions  to  our  reflected  SDE  model.  The  main 
difficulty  arises  when  we  establish  the  tightness  of  (p(°°>7)(o;,  da))^.  For  this  step,  we  focus 
on  the  particular  case  b(x)  =  —  x  (i.e.,  {Wt\t> 0  is  a  reflected  fractional  Ornstein-Uhlenbeck 
process)  where  some  explicit  computations  lead  to  a  control  of  (p(°°>7)(a;,  da))^  for  fixed 
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7  >  0.  Then  we  show  this  control  can  be  extended  to  constrained  SDE  model  whose  drift 
term  satisfies  more  general  stability  conditions.  A  joint  paper  [6]  with  Jian  Song  (University 
of  Hong  Kong)  based  on  this  topic  is  currently  under  preparation. 

3  Parameter  Estimation  for  Constrained  SDEs  Driven  by  FBMs 

3.1  Sequential  parameter  estimation  for  reflected  fractional  Ornstein-Uhlenbeck 
processes 

To  illustrate  our  problems,  consider  the  following  one-dimensional  reflected  fractional  Ornstein- 
Uhlenbeck  (RFOU)  process  on  the  positive  real  line: 

dWt  =  —r'/Wfdt  +  adB f  +  dLt , 

Wt>  0  for  all  t  >  0, 

where  7  G  (0,  00),  o  G  (0,oo)  and  BH  =  (B^)t>0  is  a  one-dimensional  fractional  Brownian 
motion  with  the  Hurst  parameter  H  G  [1/2, 1).  For  example,  the  above  reflected  fractional 
Ornstein-Uhlenbeck  process  arises  as  the  key  approximating  process  for  queueing  systems 
with  reneging  or  balking  customers  (with  long  range  dependent  arrival  process).  The  drift 
parameter  7  carries  the  physical  meaning  of  customers’  reneging  (or,  balking)  rate  from  the 
system.  In  a  recent  work  of  PI  [5],  the  goal  is 


to  estimate  the  unknown  drift  parameter  7  G  (0,  00)  based  on  observations  of  the  state 
process  {Wt}t>o- 


When  the  driving  noise  is  a  standard  Brownian  motion  (i.e. ,  H  =  1/2),  the  sequential 
estimation  plan  for  the  reflected  Ornstein-Uhlenbeck  process  has  been  studied  in  the  Pi’s 
recent  paper  [3].  More  precisely,  the  process  {Wt}  is  assumed  to  be  observed  until  the 
observed  Fisher  information  of  the  process  exceeds  a  predetermined  level  of  precision  h,  i.e., 
we  observe  {Wt}  over  the  random  time  interval  [0,  r(h)\  where  the  stopping  time  t(K)  is 
defined  as 

r(h)  :=  inf  >  0  :  J  W^ds  >  h  j>  ,  0  <  h  <  00,  (1) 

and  the  -measurable  function  7 T(h) 


lr(h)  '■ 


1 

h 


r(h) 


WAW , 


(2) 


is  a  sequential  estimator  of  the  parameter  7. 

Reccnlty,  we  obtain  the  sequential  estimation  plan  for  the  reflected  fractional  Ornstein- 
Uhlenbeck  process  ( H  G  (0, 1))  and  verify  that  the  proposed  plan  is  significantly  helpful 
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both  in  asymptotic  and  non-asymptotic  short  time  observation;  see  [5].  The  main  technical 
approach  is  based  on  the  so-called  Molchan  martingale,  that  is,  letting 

MtH  :=  [  cHs^H(t  -  sf*-H I{0<s<t}dB* 

Jo 

with  a  suitable  constant  ch  >  0,  {MfH }  is  a  Gaussian  martingale  with  square  bracket  (M)f  = 
t2~2H .  The  stopping  time  defined  by,  for  any  h  >  0, 

rH(h )  :=  inf {t  >  0  :  [  >  h}  (3) 

Jo 

plays  the  similar  role  as  in  the  standard  Brownian  motion  case  (1).  Here,  the  process  x  is  a 
suitably  constructed  functional  of  the  observable  process  W  so  that  one  can  take  advantage 
of  the  square  integrable  martingale  property  of  the  process  {/*  XsdM^  :  t  >  0}. 

With  the  stopping  time  (3),  we  construct  a  sequential  estimator  analogous  to  (2)  and 
prove  that  the  sequential  estimation  plan  (rH(h),  7^^)  has  the  following  desirable  proper¬ 
ties: 


(a)  it  is  unbiased;  (b)  the  estimation  plan  is  closed,  i.e.,  the  time  of  the  observation 
th (h)  is  finite  with  probability  1;  (c)  its  mean  squared  error  is  a  constant  that  does 
not  depend  on  the  parameter  7  to  be  estimated. 

Such  results  would  be  of  ample  use  in  applications  to  several  areas  such  as  engineering, 
financial  and  biological  modeling  where  unknown  parameter  estimation  is  based  on  relatively 
shorter  time  observation,  which  commonly  arises  in  practical  situations. 

A  paper  [5]  based  on  this  project  is  submitted  and  currently  under  review  for  Stochastic 
Processes  and  Their  Applications. 

3.2  Parameter  estimation  for  ROU  based  on  discrete  observations 

For  a  fixed  time  interval  h  >  0,  assume  that  the  discrete  time  process  {W^h  '■  1  <  k  <  n} 
are  observable  to  the  system  manager.  In  the  Pi’s  recent  work  [2],  the  following  simple 
formula  for  the  stationary  moment  of  the  reflected  Ornstein-Uhlenbeck  process  (H  =  1/2) 
is  obtained:  For  all  a  >  —1,  a  ^  0, 

E|w(co)|«  =  r  (T±b) 

where  T(-)  denotes  the  Gamma  function.  Hence,  a  natural  method  of  moment  estimator  of 
7  is  given  by 

(r((q  +  l)/2))2/°y  /  1 

7”'”'  ViEL.IwWI 
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for  an  arbitrary  h  >  0  and  any  a  >  -1,  a  ^  0.  Strong  consistency  and  asymptotic  normality 
of  the  method  of  moment  estimator  are  established  in  [2], 

The  two  main  ingredients  in  the  asymptotic  analysis  are  ergodic  properties  of  the  con¬ 
strained  stochastic  process  {W{t)}t> o  and  the  method  of  moments  estimator  constructed 
from  the  stationary  distribution  of  W(oo). 

A  paper  [2]  based  on  this  project  is  submitted  and  currently  under  review  for  Statistical 
Inference  for  Stochastic  Processes. 

3.3  Drift  parameter  estimation  for  reflected  FBM  process  based  on  the  local 
time  process 

In  the  Pi’s  recent  work  [1],  we  consider  the  one-dimensional  reflected  FBM  process  model 
with  constant  drift  parameter  b.  Suppose  the  observation  process  is  given  by  a  local  time 
process  {L(t)}t>o  (instead  of  the  state  process  {kF(t)}t>0),  which  can  be  interpreted  as  the 
total  time  the  queue  has  been  empty  up  to  time  t.  We  obtain  an  estimator  for  the  drift 
parameter  b  based  on  the  observed  cumulative  local  time  process  L,  and  establish  its  strong 
consistency  and  asymptotic  normality.  More  precisely,  for  0  <  H  <  1,  we  define  an  estimator 
bx  =  —?pL(T).  Then,  the  strong  consistency  follows: 

br  — »  b  as  T  — »  oo  a.s. 

Also,  the  asymptotic  normality  result  for  the  estimator  br  can  be  shown,  for  1/2  <  H  <  1, 

Tl~H{bT-  b)/a=>  N(0,1)  as  oo. 

The  proof  is  based  on  explicit  representations  of  the  state  and  local  time  processes  and  the 
asymptotic  results  on  the  tail  distribution  of  the  transient  queue  driven  by  FBM. 

A  paper  [1]  based  on  this  project  has  reccnlty  published  in  Journal  of  Applied  Probability. 
As  a  separate  (not  in  the  network  context)  but  remotely  related  parameter  estimation  project 
described  in  [7],  we  have  proposed  a  penalized  importance  sampling  estimator  to  obtain 
computationally  efficient  numerical  schemes  in  estimating  parameters  for  a  widely  applicable 
class  of  multidimensional  SDE  models. 

4  Non-Mar kovian  state-dependent  networks  in  critical  loading 

Queueing  systems  with  arrival  and  (or)  service  rates  depending  on  the  system’s  state  arise 
in  various  application  areas  which  include  manufacturing,  storage,  service  engineering,  and 
communication  and  computer  networks.  Longer  queues  may  lead  to  customers  being  dis¬ 
couraged  to  join  the  queue,  or  to  faster  processing,  e.g.,  when  human  servers  are  involved. 
State-dependent  features  are  present  in  congestion  control  protocols  in  communication  net¬ 
works,  such  as  TCP.  In  a  joint  paper  with  Puhalskii  [4],  we  consider  an  open  network  of  single 
server  queues  and  establish  limit  theorems  for  non-Markovian  state-dependent  networks  in 
critical  loading. 
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A  distinguishing  feature  of  the  model  is  non-Markovian  state  dependence.  More  specifi¬ 
cally,  the  number  of  customers  at  station  i,  where  i  —  1,2, ,  K,  is  governed  by  the  following 
equations 

Qi(t)  =Qi{  0)  +  Ai(t)  +  Bi(t)  —  Di(t)7 
A(t)  (J\i{Q{s))ds)  , 

JL  (4) 

3= 1 

A (t)  =NqD  ^  pi(g(s))l{Qi(s)>0}ds^  , 

where  Q(s)  =  (Qi(s), . . . ,  Qk(s ))  denotes  the  vector  of  the  queue  lengths  at  the  stations  at 
time  s . 

The  quantities  N^(t)  and  NP(t)  represent  the  number  of  exogeneous  arrivals  and  the 
maximal  number  of  customers  that  can  be  served,  respectively,  at  station  i  by  time  t  under 
“nominal”  conditions,  denotes  the  number  of  customers  routed  from  station  j  to 

station  i  out  of  the  first  m  customers  served  at  station  j  ,  and  A i(Q(t))  and  Hi(Q(t ))  represent 
instantaneous  exogenous  arrival  and  service  rates,  respectively,  for  station  i  at  time  t  given 
the  queue  length  vector  Q(t) .  Generalised  Jackson  networks  is  a  special  case  of  (4)  where 
the  Nf  and  Nf  are  renewal  processes,  the  Ty  are  Bernoulli  processes,  and  A.,  ( • )  =  /i,  ( • )  =  1 . 

We  obtain  limit  theorems  in  critical  loading  for  the  queue  length  processes  akin  to  clas¬ 
sical  diffusion  approximation  results  available  for  generalised  Jackson  networks.  The  limit 
stochastic  process  is  a  continuous-path  reflected  process  (generally,  non-Markovian)  on  the 
nonnegative  orthant.  The  constrained  stochastic  differential  equation  driven  by  FBM  is  a 
special  example  of  such  limit  stochastic  processes.  We  give  an  application  to  generalised 
Jackson  networks  with  state-dependent  rates. 

The  main  contribution  of  our  work  [4]  is  incorporating  general  arrival  and  service  pro¬ 
cesses.  This  is  achieved  by  applying  an  approach  different  from  the  ones  used  by  other 
related  works.  In  our  approach  we,  in  a  certain  sense,  return  to  the  basics  and  employ  ideas 
which  have  proved  their  worth  in  the  set-up  of  generalised  Jackson  networks.  We  show  that 
continuity  considerations  may  produce  stronger  conclusions  at  less  complexity.  Our  main 
result  states  that  if  the  network  primitives  satisfy  certain  limit  theorems  with  continuous- 
path  limits,  then  the  multidimensional  queue-length  processes,  when  suitably  scaled  and 
normalised,  converge  to  a  reflected  continuous-path  process  on  the  nonnegative  orthant. 

A  paper  [4]  based  on  this  project  is  submitted  and  currently  under  review  for  Stochastic 
Models. 
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